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In order to understand the conditions which lead a highly magnetized, relativistic plasma to 
become unstable, and in such cases how the plasma evolves, we study a prototypical class of magne¬ 
tostatic equilibria where the magnetic field satisfies V x B = ctB, where a is spatially uniform, on 
a periodic domain. Using numerical solutions we show that generic examples of such equilibria are 
unstable to ideal modes (including incompressible ones) which are marked by exponential growth in 
the linear phase. We characterize the unstable mode, showing how it can be understood in terms of 
merging magnetic and current structures, and explicitly demonstrate its instability using the energy 
principle. Following the nonlinear evolution of these solutions, we find that they rapidly develop 
regions with relativistic velocities and electric fields of comparable magnitude to the magnetic field, 
liberating magnetic energy on dynamical timescales and eventually settling into a configuration 
with the largest allowable wavelength. These properties make such solutions a promising setting for 
exploring the mechanisms behind extreme cosmic sources of gamma rays. 


Introduction .—Magnetic stability is a fundamental 
question in a range of fields from laboratory plasma 
physics, where it influences the viability of fusion de¬ 
vices [1]; to space physics, where it controls the struc¬ 
ture of magnetic fields within stars and planets [2]. In 
high-energy astrophysics, the spontaneous release of en¬ 
ergy associated with transitions between magnetic equi¬ 
librium states is of particular importance to understand¬ 
ing the dramatic gamma-ray activities from pulsar wind 
nebulae [3, 4], magnetars [5-8], relativistic jets associated 
with active galactic nuclei [9-12], and gamma-ray bursts. 
These diverse sources exhibit powerful gamma-ray flares 
on timescales short compared with their light-crossing 
times [7, 11, 12], and seem to require that electrons and 
positrons be accelerated throughout extended regions, to 
energies as high as several PeV [13, 14]. The most dra¬ 
matic variations are likely produced in the relativistic 
electromagnetic outflows away from the central engine 
(neutron stars or black holes), and a mechanism is press- 
ingly needed to explain the rapid, volumetric conversion 
of magnetic energy into high energy particles and radia¬ 
tion. Here, we consider whether such a process may be 
triggered by magnetic instability in the outflow. These 
outflows may initially accelerate, so that they cannot be 
crossed by hydromagnetic waves in an outflow timescale. 
However, they will eventually be decelerated when their 
momentum flux decreases to that of the external medium, 
bringing disconnected regions back into causal contact 
where they are likely to be unstable ^ . 

To understand the conditions under which a plasma 
becomes unstable, and to follow its subsequent nonlinear 
evolution in an idealized setting, we focus on a model 
class of force-free equilibria, which we find evolves in a 


^ A cosmological analogy would be when perturbations from the 
epoch of inflation are believed to have “re-entered” the horizon 
and exhibited gravitational instability. 


manner that is both surprising on formal grounds, and 
highly suggestive of the behavior of the most dramatic 
cosmic sources. Force-free solutions, where the Lorentz 
force vanishes, are an excellent approximation for highly 
conducting and strongly magnetized plasmas, where the 
plasma inertia and pressure is sub-dominant to the mag¬ 
netic field, and have been used extensively across differ¬ 
ent fields. A particularly important class of force-free 
equilibria that are conjectured to arise naturally from 
magnetic relaxation are the so called Taylor states, which 
satisfy the Beltrami property: V x B = aB where a 
is a constant [15]. These solutions have played an im¬ 
portant role not only in laboratory plasma physics [16], 
but also in solar physics [17-19], astrophysics [20], and 
beyond [21]. In this work we focus on space-periodic 
equilibria as a simple, computationally tractable setting 
free of the effect of confining boundaries (as in extended 
outflows). Though there is a rich literature studying such 
solutions [15, 22-25], important facts regarding their sta¬ 
bility have not been appreciated. Focusing on a prototyp¬ 
ical example, the “ABC” solutions [26] (defined below), 
in [24] it was claimed that such solutions are stable to in¬ 
compressible perturbations (see also [25]). Here we show 
that, in fact, generic periodic Beltrami magnetic fields 
are linearly unstable, including to incompressible defor¬ 
mations. The only exceptions we find are special cases 
lacking magnetic curvature, and those in the fundamen¬ 
tal mode or ground state, having the lowest magnetic 
energy compatible with conservation of magnetic helic- 
ity Hm = f a ■ 'BdV (where A is the magnetic vector 
potential). The instability we find is ideal, in contrast 
to previous studies of dissipative effects [27]. We find 
that in the nonlinear evolution, magnetic energy is in¬ 
deed liberated rapidly, giving rise to relativistic veloci¬ 
ties and electric fields of comparable magnitude to the 
magnetic fields on dynamical timescales, and eventually 
allowing the system to relax to its ground state. These 
solutions are therefore a simple, but promising setting 
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to explore the mechanisms underlying extreme cosmic 
sources of gamma rays. 

In what follows, we present simulation results show¬ 
ing the linear-regime instability of a range of magneto¬ 
static equilibria, and then illustrate the properties of the 
dominant unstable mode in some example cases, inde¬ 
pendently confirming the growth rate using the energy 
principle. We then compare the results found using var¬ 
ious degrees of magnetization, discuss the nonlinear evo¬ 
lution of the instability, and conclude. We use units with 
c = 1 throughout. 

Methodology. —The equilibrium magnetic fields we 
study are exemplified by the three-parameter “ABC” 
field [26, 28] given by 

= (^3 cos az — B 2 sin ay^ ( 1 ) 

Bi cos ax — Bs sin az, B 2 cos ay — Bi sin ax ). 

We use some particular examples of this equilibrium solu¬ 
tion for illustrative purposes, but also consider the more 
general class of Beltrami fields [29] B = a^ + V x ^ 
where the potential field ^ is any solenoidal vector field 
satisfying the vector Helmholtz equation = 0, 

so that ^ comprises only the Fourier harmonics whose 
wave-vector k has magnitude a. These more general con¬ 
figurations are constructed by choosing random vector 
amplitudes for the admissible harmonics. Our compu¬ 
tational domain is the periodic cube of length L = 27r 
(though we restore L in some places for clarity). 

We simulate a perfectly conducting, magnetized 
fluid and consider cases with different finite values 
of the volume-averaged magnetization parameter a := 
{B‘^/4:7Tph) where ph is the fluid enthalpy (treated us¬ 
ing the ideal relativistic magnetohydro dynamic equations 
and the code in [30]), as well as the limiting case of a com¬ 
pletely magnetically dominated plasma, cr = oc (treated 
by force free electrodynamics [31-33]). See supplemental 
material below for details. 

Instability in the Linear Regime. —For this class of 
magnetic equilibria, we find generic solutions with a^ > 1 
to be unstable to linear ideal perturbations that are 
characterized by exponential growth of the electric field. 
Fig. 1 illustrates this for a case with a^ = 11. (Here we 
present results from a = 00 simulations, and in a later 
section we compare these to finite magnetization cases.) 
The magnitude of the growing solution is proportional to 
the initial perturbation^, consistent with a linear insta¬ 
bility. The growth rate of, e.g., the electric field energy 
is converging to 7 ~ A.Oa/L with increased resolution, 
evidence that the instability is not due to numerical/non¬ 
ideal effects. 


^ Here the initial perturbation we use is an electric field with = 
Eq cos( 27 Ty/L) and where the other components are given by 
cyclic permutations of {x,y, z}. From this we subtract out any 
component parallel to B. 



FIG. 1. Results from simulations with = 11. Top: 
The growth of the electric field energy He, normalized by 
its initial value, for three different values of the initial per¬ 
turbation. Middle: The growth in He for three different res¬ 
olutions, along with an exponential fit. The difference be¬ 
tween the best fit exponent for the high resolution, and the 
Richardson extrapolated value using all three resolutions, is 
~ 0.1%, the extrapolation being consistent with between first 
and second-order convergence. The bottom panel illustrates 
the conservation of total energy U for three different resolu¬ 
tions. Though initially higher-order when the equilibrium- 
solution truncation error dominates, the convergence even¬ 
tually drops to first-order, presumably because (as discussed 
below) the unstable solution has non-smooth features. Con¬ 
servation of magnetic helicity is similar. 

Other equilibrium solutions exhibit similar exponen¬ 
tially growing solutions, as shown for some example cases 
in Fig. 2. This holds for wavelengths larger than the fun¬ 
damental mode for the domain, a^ = 1, which is known 
to be stable [24, 25]. The growth rate of the instability 
is also roughly proportional to a, though there is depen¬ 
dence on the particular realization used. 

The Dominant Unstable Mode .—In order to illustrate 
the nature of the instability, we focus on a simple type of 
a = 2 equilibrium solution given by Eq. 1 . We illustrate 
this solution for three different choices of coefficients in 
Fig. 3. As discussed in [28] for the mathematically equiv¬ 
alent Euler flow, these solutions have a rich structure. 
The {Bi, B 2 , Bs) = (1,1,0) case consists of “vortices:” 
regions of helical field (and current) lines circling a cen¬ 
tral axis. The {Bi, B 2 , B^) = (1,1/2,0) case has vor¬ 
tices as well as “shear layers:” wavy field lines that begin 
and end on opposite sides of the domain. In addition to 
these two cases with z-translational symmetry, we also 
show a more generic case where all three coefficients are 
nonzero (and that like the second case, has no places 
where = 0 ). 

In Eig. 3 we also show the corresponding velocity field 
V = E X B^/|B^p (which will be proportional to the dis- 
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FIG. 2. The growth in electric field energy for various values 
of and three different realizations for = 11. Time has 
been shifted so that all the curves have an abscissa of 0 at the 
ordinate value of 10“®, and the time axis has been scaled by 
a which gives the different examples roughly the same slope. 
The = 1 simulation does not exhibit exponential growth, 
and has been scaled up by an overall factor. 


placement ^ for an eigenmode) characterizing the dom¬ 
inant instability arising in each case. This is calculated 
from a numerical snapshot after the instability has grown 
by roughly 10 orders of magnitude — seeded in this case 
just by truncation error — but is still in the linear regime 
(|E| ^ The velocity field acts to bring to¬ 

gether vortices, or current channels, circulating in the 
same direction in order to move towards a larger wave¬ 
length, lower magnetic energy configuration. The nonlin¬ 
ear evolution is characterized by the merging of magnetic 
vortices and can thus be related to the coalescence insta¬ 
bility of magnetic islands [38]. We can also see that the 
velocity field appears to have non-smooth features, rem¬ 
iniscent of spontaneous current sheets [39], that occur at 
the separatrices dividing the vortices and shear layers. 
Though the generic case lacks ^-translational symmetry, 
it appears qualitatively similar to the second case. 

From the (5i, ^ 2 , Bs) = (1,1, 0) to the (1,1/2, 0) case, 
the growth rate of the instability decreases by a factor of 
~ 1.9 with the addition of the shear layers in the equilib¬ 
rium solution. In fact, the growth rate decreases mono- 
tonically with B 2 , and as ^2 0 and the vortices shrink 

to zero volume, the growth rate of the instability also 
goes to zero. In fact, it can be shown (see supplemen¬ 
tal material below) that the single mode solutions are all 
stable, including those at short wavelengths. 

For the generic case where the three coefficients are 
nonzero, we directly confirm that the instability is linear 
and ideal by using the energy principle [40] . This implies 
that an equilibrium solution satisfying the Beltrami 


property is unstable to a displacement ^ if 
2 _ /[|(5B|2-a(exBq.(5B)]dy 

where (5B := V x x B^), and that the instability should 
grow at least as fast as jcjl. Computing this quantity 
using finite differences, with the numerical velocity field, 
gives a value of \uj\L ^ 2.4, near the growth rate of 2.6 
measured for the electric field in the simulation. 

Although the unstable displacement from the simu¬ 
lations has both non-smooth and compressible (/ |V • 
v\dV ^ 0.15 /|V X v\dV) features, these are not nec¬ 
essary for instability. For example, by applying a low- 
pass filter to the Fourier harmonics of ^ one can con¬ 
struct an ideal perturbation that has no power at scales 
|k| > 29, and yet gives a value of \uj\L ^ 2.3 upon eval¬ 
uating Eq. ( 2 ) algebraically in /c-space. Furthermore, we 
have also explicitly constructed (see supplemental ma¬ 
terial below) analytical examples of smooth and incom¬ 
pressible displacements that destabilize particular Bel¬ 
trami solutions, which are counterexamples to the claims 
in [24]. 

Magnetization comparison and nonlinear evolution .— 
The same qualitative behavior is observed when the 
plasma is evolved with different magnetization param¬ 
eters. Fig. 4 shows the evolution of the kinetic and mag¬ 
netic energy Ub for different parameters a values, in¬ 
cluding the limiting case of a = oc. The lower inset of 
Fig. 4 shows that the growth rate of the unstable mode 
increases monotonically with increasing a and is roughly 
proportional to the Alfven speed va = \/(Jjiy + cr). 

In both finite and infinite magnetization cases, expo¬ 
nential growth of the unstable displacement persists until 
its velocity |v| approaches the Alfven speed va (near the 
speed of light for large a) and higher-order evolutionary 
terms dominate. Following the turbulent state, the sys¬ 
tem settles into a lower energy equilibrium. Somewhat 
surprisingly, we do not observe direct evolution into the 
lowest energy {a = 1) state for all cases. For example, 
the = 11 state in Fig. 4 first transitions into a config¬ 
uration with ^ 97% of its spectral energy in modes with 
= 3, where it remains for about ten Alfven crossing 
times, before making a second transition into the ground 
state, where 99% of its energy in k^ = 1 modes. The 
lifetime of the intermediate state may be related to the 
fact that ^ 88 % of the energy is in a single k^ = 3 mode. 

During the nonlinear evolution, regions develop where 
|E| is comparable, and even exceeds |B|. Since the hyper- 
bolicity of the equations breaks down for cr = 00 when 
this happens, to evolve further we handle such regions 
with an ad-hoc prescription where we simply reduce the 
electric-field magnitude to equal that of the magnetic 
field, leading to a reduction of energy. The finite a cases 
do not suffer from this problem, and our scheme explicitly 
conserves total energy, but permits conversion of mag¬ 
netic or kinetic energy into internal energy, especially at 
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FIG. 3. Streamlines of a magnetic field equilibrium solution given by Eq. 1 with ol — 2 and various coefficients (top), and the 
corresponding velocity field v = E x of the unstable mode arising from the simulations (bottom) in the z = 0 plane. 

The equilibrium solutions, from left to right, correspond to (Bi,B 2 ,B 3 ) = (1,1,0), (1,1/2, 0), and (—0.814,0.533,0.232), 
respectively. The color indicates the perpendicular vector component with red and blue representing, respectively, out of the 
page and into the page. The thickness of the streamline is proportional to the vector magnitude. The black lines indicate the 
location of the separatrices in the equilibrium solutions. 


shocks or places where the magnetic field is nearly dis¬ 
continuous. Encouragingly, we still find consistency be¬ 
tween these different (and somewhat arbitrary) types of 
energy dissipation in the non-linear regime. For example, 
as shown in Fig. 4, we find the same energy levels associ¬ 
ated with the intermediate and final magnetic equilibria. 
This is consistent with conservation of magnetic helicity. 
Since the Beltrami fields have B = a A, their helicity is 
2UbI(^^ and the ratio of magnetic energy in the ai and 
af equilibria is simply afjai. Accordingly, we do not ex¬ 
pect the dissipation mechanism to have much influence 
on the energy of the final state if helicity is preserved. 
(For the simulations shown in Fig. 4, Hm is constant to 
^ 0.1%.) But understanding the details of the energy 
dissipation will require better physical modeling. 

Conclusions .—We studied periodic Beltrami magnetic 
fields and found that they were unstable to ideal modes. 
Though we focused on the relativistic case, since the in¬ 
stability is linear in velocity, it also applies to the nonrel- 
ativistic setting. This contrasts with [24] — where it was 
concluded that such solutions are linearly stable against 
incompressible perturbations — and suggests that the re¬ 
laxation of a complex magnetic field will not terminate 
at a small wavelength equilibrium, but instead undergo 
a so-called inverse helicity cascade where magnetic en¬ 
ergy goes to the largest available scale [41-44]. There 



FIG. 4. A comparison of the decay of an = 11 equilibrium 
in simulations with different values of magnetization parame¬ 
ter a. Shown is the magnetic energy (top) and kinetic/electric 
field energy (bottom). The horizontal dashed lines in the top 
panel indicate the magnetic energy of = 3 and = 1 
states with the same helicity. The bottom inset shows the 
linear growth rate 7 measured for runs having different mag¬ 
netization parameters, along with the Alfven speed (dashed 
line) for comparison. 
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are known examples of unstable cylindrically-symmetric 
Beltrami solutions [45], and studying a broader class of 
geometries would be interesting follow-up work. 

For highly magnetized, relativistic plasma, the insta¬ 
bility gives rise to regions where the electric field mag¬ 
nitude is comparable to the magnetic field on dynamical 
timescales. In extreme cosmic sources of gamma rays, 
where such configurations may be relevant, these would 
be likely sites of particle acceleration and photon emis¬ 
sion. Understanding the details of this, including the role 
of magnetic reconnection [46-48] and turbulence [49] in 
ultimately dissipating energy, and determining the na¬ 
ture of the acceleration mechanism [50-52], will require 
kinetic simulations incorporating radiative losses, some¬ 
thing we plan for future work. 
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Supplemental material: Details of numerical methods 


In the main text we assume a perfectly conducting medium and perform simulations of cases having various degrees 
of magnetic dominance. The limiting case of a completely magnetically dominated plasma is treated by FFE [31-33]. 
The evolution equations of FFE are just the Maxwell equations with a prescription for the current derived from the 
assumption that the Lorentz force vanishes and can be written as [33, 34] (we use Heaviside-Lorentz units and set 
c = 1 throughout): 


dB _ _ 

j = T [(B . (V X B) - E ■ (V X E)) B + (V ■ E)E X B]. 


(3a) 

(3b) 

(3c) 


We numerically solve these using fourth-order finite difference stencils and Runge-Kutta time stepping. We use 
hyperbolic divergence cleaning to exponentially damp violations of the V • B = 0 constraint as in [35] . The E • B = 0 
constraint is explicitly enforced by redefining E ^ E — B(E • B)/H^ at every coarse time step. We apply standard 
sixth-order Kreiss-Oliger [36] numerical dissipation to all the hyperbolic variables to suppress high frequency numerical 
error. The FFE code is parallelized using the PAMR/AMRD software library 

We also solve the ideal relativistic magnetohydro dynamic (RMHD) equations (with a F = 4/3 equation of state) for 
different values of the volume-averaged magnetization parameter a := {B‘^ jAnph) (where ph is the fluid enthalpy). In 
addition to the specified magnetic field, we use a restmass density and pressure that are initially equal and uniform. We 
use a second-order, constrained-transport, finite-volume scheme that explicitly conserves mass, energy, momentum, 
and magnetic flux. Full details of the code are described in [30]. 


Supplemental material: Unstable modes with analytical methods 


In this sections, we verify the existence of unstable ideal modes for periodic Beltrami solutions using analytical 
methods. Introducing v = E x B/H^ and E = — v x B, the evolution equations (3) can be rewritten in terms of the 
following equations [37] 


dt 


V X (v X B), 


m 


[B\) = (V X B) X B + (V X E) X E + (V • E)E 


(4a) 

(4b) 


^ http://laplace.physics.ubc.ca/Group/Software.html 
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Now consider small perturbation on a static equilibrium state which satisfies V x Bq = aBo with a being a constant. 
Let B = Bq + Bi, and define the displacement field ^ such that to first order d^/dt = v. Suppose the perturbation 
can be decomposed into normal modes oc The linearized version of Equation (4) now becomes 

Bi = V X (|*x Bo), (5a) 

—= (V X Bi — ofBi) X Bq 

= (V X Bo) X (V X (f X Bo)) + [V X (V X (f x Bo))] x Bo (5b) 

In the equation for ^ the differential operator K is self-adjoint under periodic boundary conditions [40], so we can 
define a potential energy 

y = -\ y'd3^{[Vx{fxB)]2-a(exB).[Vx(exB)]} (6) 


As a result 


This allows a variational approach to the stability problem: we can use trial functions to get an upper limit on the 
lowest An equilibrium state is unstable as long as we can find one trial function that renders V negative. In the 
following we show a few examples. 


ID equilibria 


One-dimensional force-free equilibria in periodic box can be written as Bq = {cos sin 2 p{z), 0), where ^{z) = az 
for linear force-free fields. We can write the normal mode perturbation as ^ component 

of ^ parallel to the background Bq does not have physical significance, we can impose the requirement ^ _L Bq. Then 
from Equation (5), we get 


ix{z) 

d 


i sin 7 /; {ky cos ip — kx sin 2 p) , 


kl ^ k‘^ — uo‘^ 




— {kx COS Ip -h ky sin?/^) ^ dz 


^(cj^ - {kx cosp; -h ky smp>) ^) -h (cj^ - - ky) = 0 


(8a) 

(8b) 


Eigensolutions with uj‘^ > k^ ky exist and they correspond to oscillating normal modes. However, we do not have 
normal modes with < 0: otherwise — {kx cos 2 p ky sin p>)‘^ <0 and the solution would be essentially exponential 
so it cannot satisfy the boundary condition. As a result, the ID equilibria with periodic boundary conditions are 
stable to ideal modes. 


2D and 3D equilibria 


In this case it’s not realistic to solve the normal mode equation so we make use of the variational principle. In periodic 
box it’s convenient to use Eourier basis to construct our parameterized trial functions: ^ where 

k_^ = —k^, = C* fo ensure reality. Bq can also be decomposed into Eourier components: Bq = 

where a_n = —a^, B_n = B*, • B^ = 0, ja^l = \a\ and ia^ x B^ = o^B^. Then the integrals for calculating the 
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potential energy in Equation (7) only involve algebraic manipulations (cf [24], with modifications): 

n2Tv /*27r /*27r 


V = 


\f n 


Bi • [Bi — ^ X (V X Bo)] dxdydz = (Bi • [Bi — ^ x (V x Bq)]) 


— ^ ^ ^ ^ ^ ^n') ^ (Cm' ^ ^n')] ’ [(^m "b Q-n) ^ (Cm B^) X (aT^ X B^)] 


K=k^+a^ n n' 
=k^+a^ 


E U 


K=k^+a, 

1 


E[(k m B„)^ m O^m ' Cm)^^ 


^ ^ (Cm ' Bn)^r) 


^ ^ (Cm ' ^n)Br 


(9) 


^^[(km ’ B 72 )Cm (km ’ Cm)Bn (Cm ’ ^n)^n (Cm ' 


For certain equilibria this expression is not positive definite and can turn out to be negative using appropriate trial 
functions. 

In the following we focus on a particular case, the ABC field 

B = 5i(0, cosax, — sin ax) + B2{— sin ay, 0, cosa^) + 53(cosa2:, — sinaz, 0), (10) 


which has been studied in the main text using numerical simulations. 

The first case is 2D with Bi = B 2 = 1, B^, = 0,a = 2 > 1. We include several of the longest wavelength 
perturbations in the trial function: k= (1,1,0), (1,—1,0), (1,3,0), (1,—3,0), (3,1,0), (3,—1,0), (3,3,0), (3,—3,0) 
plus their negative companions, then decide the coefficients by minimizing the right hand side of Equation (7). The 
k’s are chosen based on the observation that in the numerical simulation, the dominant unstable mode is 2D and has 
an electric field E oc ^ x B mainly comprised ofk= (±1,±1,0) components. We find that the minimization does give 
negative < —0.04, meaning a growth rate of 7 > 0 . 2 . As a comparison, the light crossing time scale r = 27r so 

yr > 1.26. This growth rate is less than what has been found from numerical simulations (yr ^ 5.5 for this particular 
case) which is to be expected since we can only get a lower limit on the growth rate from a variational approach. The 
Fourier components of this trial function are listed in Table I, and Figure 5 shows the stream plot of the perturbation. 
We find that including higher k’s (but only k = (±(2s + l),±(2t + l), 0 ), 5 ,t G Z are relevant) allows us to get 
lower minimum (i.e. larger growth rate). For example, when the Fourier modes for the perturbation ^ include 
the following: k = (1,1,0), (1,-1,0), (1,3,0), (l,-3,0), (3,1,0), (3,-1,0), (3,3,0), (3,-3,0), (1,5,0), (l,-5,0), 
(5,1,0), (5,—1,0), (3,5,0), (3,—5,0), (5,3,0), (5,—3,0), (5,5,0), (5,—5,0), we get the growth rate y > 0.38 and 
yr > 2.38. The corresponding trial function is plotted in Figure 6 . Thus, the short wavelength perturbations are 
essential for the instability. Comparing these analytical trial functions with the dominant unstable mode discovered 
in numerical simulations, we find that the former is consistent with being truncated version of the latter. 


TABLE I. Fourier components of the unstable trial function (com¬ 
pressible) for the equilibrium Bi — B 2 — B^ — a — 2 

{1,1, 0} {0.0795, 0.0688, 0.0641} {0.1367, 0.0116, 0.0502} 

{1,-1,0} {-0.0623,-0.0569,0.1211} {-0.1611,-0.1557,0.0766} 
{1, 3, 0} {-0.0646, 0.0145, 0.0474} {0.0342, 0.0145, 0.0515} 

{1, -3, 0} {-0.0400, 0.0045, 0.1013} {-0.0972, 0.0045, -0.0441} 
{3,1, 0} {-0.0145, 0.0700, -0.0420} {0.0145, 0.0289, 0.0568} 

{3, -1, 0} {0.0045, 0.1083, 0.0470} {-0.0045, -0.0511, 0.1042} 

{3, 3, 0} {0.0048, -0.0048, 0.0008} {-0.0048, 0.0048, 0.0008} 

{3, -3, 0} {-0.0154, -0.0154, 0.0027} {0.0154, 0.0154, 0.0027} 

< —0.04, T — 27r, UT > 1.26 

^ The coefficients for the —km terms are just the complex conju¬ 
gate of the listed km coefficients. 


We also find that the instability still exist when incompressibility is imposed, i.e. V • ^ = 0. Table II lists the 
Fourier components for the incompressible trial function that gives negative potential energy and Figure 7 shows the 
corresponding stream plot. 
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FIG. 5. (a) The equilibrium ABC field (10) with = B 2 = 1 , B 3 = 0 , a = 2 > 1 . The stream lines indicate the field 
components in the x — y plane and color indicates the component perpendicular to the plane (same below for other vector 
fields), (b) Trial perturbation ^ that renders negative potential energy (instability) for the equilibrium. We only plot the 
components perpendicular to the equilibrium magnetic field. The perturbation is compressible in this case, (c) Perturbation 
magnetic field Bi resulted from the perturbation in (b). 



FIG. 6 . (a) Similar to Figure 5 (b), this is the displacement field derived from a trial function involving more Fourier modes, 
(b) Perturbation magnetic field Bi resulted from the perturbation in (a). 


Other more general 2 D cases with Bi ^ B 2 and B 1 B 2 7 ^ 0 are found to be unstable as well in numerical simulations, 
with growth rate decreasing from maximum to 0 as, say, B 2 IB 1 goes from 1 to 0 . As an illustration we apply our 
variational principle to the case Bi = 1, B 2 = 1/2, Bs = 0, using trial functions comprised of Fourier modes 
k = (1,1, 0), (1, —1, 0), (1, 3, 0), (1, —3, 0), (3,1, 0), (3, —1, 0), (3, 3, 0), (3, —3, 0) and their negative companions. This 
is shown in Table III Figure 8 , and the displacement field can be readily compared with Figure 3 in the main text. 
The growth rate we get from this trial function is yr > 0.8 — the lower limit has been reduced from corresponding 
Bi = B 2 = I case. 

It is also instructive to consider a 3D example. Here we take Hi = H 2 = 1, H 3 = 1/5 and a = 2. Still use Fourier 
components k = (1,1, 0), (1, —1,0), (1,3, 0), (1, —3, 0), (3,1, 0), (3, —1, 0), (3, 3, 0), (3, —3, 0) in the trial function, we 
found the unstable perturbation as shown in Table IV and Figure 9. 
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FIG. 8 . (a) The equilibrium ABC field ( 10 ) with Bi = 1, B 2 = 1 / 2 , B 3 = 0 , a = 2 > 1 . The stream lines indicate the field 
components in the x — y plane and color indicates the component perpendicular to the plane (same below for other vector fields). 
The thick black lines show the separatrices — surfaces separating different topological domains — in the equilibrium magnetic 
field, (b) Trial perturbation ^ that renders negative potential energy (instability) for the equilibrium. We only plot ^j_, the 
components perpendicular to the equilibrium magnetic field. The perturbation is compressible in this case, (c) Perturbation 
magnetic field Bi resulted from the perturbation in (b). 
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TABLE IV. Fourier components of the unstable trial function (com¬ 
pressible) for the equilibrium Bi — B 2 — B 3 = 1/5, a — 2 
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FIG. 9. (a) The equilibrium ABC field (10) with Bi — B 2 — 1, B 3 = 1/5, a — 2 > 1 . The stream lines indicate the field 
components in the z = 0 plane and color indicates the component perpendicular to the plane (same below for other vector 
fields), (b) Trial perturbation ^ that renders negative potential energy (instability) for the equilibrium. We only plot the 
components perpendicular to the equilibrium magnetic field. The perturbation is compressible in this case, (c) Perturbation 
magnetic field Bi resulted from the perturbation in (b). Both (b) and (c) are stream plots on the z = 0 plane. 
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